rds=r(0)*rr     ;m
rdsi=1./rds
GG=6.674e-11   ;N(m/kg)^2
mass=1.988e30     ;kg
dz=(r[l-1]-r[0])*rr/double(l-1)
g=GG*mass/(rr*r)^2     ;m/s^2
g0=g[0]
z = indgen(l)*dz

Rg=13732.
cp=2.5*rg
cap=rg/cp
capi=1./cap
gamma=cp/(cp-rg)

if (itac eq 1) then begin 
 ;;BOTTOM
 rho00 = 370.33
 T00 = 3.4976e6
 p00 = 1.7788e13 
 th00= T00
 g0=734.428
 kk=0
;
endif else begin
;;BASE OF CZ (l=35)
 rho00=167.74
 p00=4.8454e12
 T00=2.10357e6
 th00 = T00
 g0=528.35
endelse
;
Rg=13732.
cp=2.5*rg
cap=rg/cp
capi=1./cap
ss=g0/(cp*th00)
rh = rho00*(1.-ss*(z-z[0])/(1.+(z-z[0])*rdsi))^(capi-1)
pr=p00*(rh/rho00)^(gamma)
temp=pr/(rh*Rg)
icz=where(r ge 0.75)
cs = sqrt(gamma*Rg*mean(temp[icz]))

rho_av=mean(rh)

rh2d = fltarr(l,m)
T2d = fltarr(l,m)

for k=0, l-1 do begin
  for j=0,m-1 do begin
    rh2d[k,j]=rh[k]
;    rh2d[k,j]=1.
    T2d[k,j]=temp[k]
  endfor
endfor
